{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10, 5]\n",
    "\n",
    "import pdb\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## inspired by https://rajeshrinet.github.io/blog/2014/ising-model/\n",
    "\n",
    "### CONFIGURATION\n",
    "N           = 5 # width of lattice\n",
    "M           = 5 # height of lattice\n",
    "temperature = 0.5\n",
    "BETA        = 1 / temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialize system\n",
    "def initRandState(N, M):\n",
    "    block = np.random.choice([-1, 1], size = (N, M))\n",
    "    return block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def costForCenterState(state, i, j, n, m):\n",
    "    centerS = state[i, j]\n",
    "    neighbors = [((i + 1) % n, j), ((i - 1) % n, j),\n",
    "                 (i, (j + 1) % m), (i, (j - 1) % m)]\n",
    "    ## notice the % n because we impose periodic boundary conditions\n",
    "    ## ignore this if it doesn't make sense - it's merely a physical constraint on the system\n",
    "\n",
    "    interactionE = [state[x, y] * centerS for (x, y) in neighbors]\n",
    "    return np.sum(interactionE)\n",
    "\n",
    "\n",
    "def magnetizationForState(state):\n",
    "    return np.sum(state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mcmc steps\n",
    "def mcmcAdjust(state):\n",
    "    n = state.shape[0]\n",
    "    m = state.shape[1]\n",
    "    x, y = np.random.randint(0, n), np.random.randint(0, m)\n",
    "    centerS = state[x, y]\n",
    "    cost = costForCenterState(state, x, y, n, m)\n",
    "    if cost < 0:\n",
    "        centerS *= -1\n",
    "    elif np.random.random() < np.exp(-cost * BETA):\n",
    "        centerS *= -1\n",
    "    state[x, y] = centerS\n",
    "    return state\n",
    "    \n",
    "def runState(state, n_steps, snapsteps = None):\n",
    "    if snapsteps is None:\n",
    "        snapsteps = np.linspace(0, n_steps, num = round(n_steps / (M * N * 100)), dtype = np.int32)\n",
    "    saved_states = []\n",
    "    sp = 0\n",
    "    magnet_hist = []\n",
    "    for i in range(n_steps):\n",
    "        state = mcmcAdjust(state)\n",
    "        magnet_hist.append(magnetizationForState(state))\n",
    "        if sp < len(snapsteps) and i == snapsteps[sp]:\n",
    "            saved_states.append(np.copy(state))\n",
    "            sp += 1\n",
    "    return state, saved_states, magnet_hist\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### RUN A SIMULATION\n",
    "init_state = initRandState(N, M)\n",
    "plt.imshow(init_state)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_state, states, magnet_hist = runState(init_state, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(final_state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(magnet_hist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Modify the simulation code above to visualize magnetization over time for 100 test runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = []\n",
    "for i in range(100):\n",
    "    init_state = initRandState(N, M)\n",
    "    final_state, states, magnet_hist = runState(init_state, 1000)\n",
    "    results.append(magnet_hist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for mh in results:\n",
    "    plt.plot(mh,'r', alpha=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: generate a curve of absolute value of magnetization at step 100 for temperature for temperatures = 0.1, 0.5, 1, 2, 5, 20, 50, 100, 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = []\n",
    "temps = [0.1, 0.5, 1, 2, 5, 20, 50, 100, 500]\n",
    "\n",
    "# mcmc steps\n",
    "def mcmcAdjust(state, beta):\n",
    "    n = state.shape[0]\n",
    "    m = state.shape[1]\n",
    "    x, y = np.random.randint(0, n), np.random.randint(0, m)\n",
    "    centerS = state[x, y]\n",
    "    cost = costForCenterState(state, x, y, n, m)\n",
    "    if cost < 0:\n",
    "        centerS *= -1\n",
    "    elif np.random.random() < np.exp(-cost * beta):\n",
    "        centerS *= -1\n",
    "    state[x, y] = centerS\n",
    "    return state\n",
    "    \n",
    "def runState(state, beta, n_steps, snapsteps = None):\n",
    "    if snapsteps is None:\n",
    "        snapsteps = np.linspace(0, n_steps, num = round(n_steps / (M * N * 100)), dtype = np.int32)\n",
    "    saved_states = []\n",
    "    sp = 0\n",
    "    magnet_hist = []\n",
    "    for i in range(n_steps):\n",
    "        state = mcmcAdjust(state, beta)\n",
    "        magnet_hist.append(magnetizationForState(state))\n",
    "        if sp < len(snapsteps) and i == snapsteps[sp]:\n",
    "            saved_states.append(np.copy(state))\n",
    "            sp += 1\n",
    "    return state, saved_states, magnet_hist\n",
    "\n",
    "\n",
    "res = []\n",
    "for temp in temps:\n",
    "    temp_res = []\n",
    "    for _ in range(20):\n",
    "        init_state = initRandState(N, M)\n",
    "        final_state, states, magnet_hist = runState(init_state, 1/temp, 100)\n",
    "        temp_res.append(abs(magnet_hist[-1]))\n",
    "    res.append(np.mean(temp_res))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(temps, res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: what might an agent-based model look like when applied to this problem?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
